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Foreground removal is a challenge for 21-cm tomography of the high redshift Uni- 
verse. We use archival GMRT data (obtained for completely different astronomical 
goals) to estimate the foregrounds at a redshift ~ 1. The statistic we use is the cross 
power spectrum between two frequencies separated by Av at the angular multipole 
I, or equivalently the multi-frequency angular power spectrum Ct(Av). An earlier 
measurement of Ce(Av) using this data had revealed the presence of oscillatory pat- 
terns along Av, which turned out to be a severe impediment for foreground removal 
(Ghosh et al. 2011). Using the same data, in this paper we show that it is possible 
to considerably reduce these oscillations by suppressing the sidelobe response of the 
primary antenna elements. The suppression works best at the angular multipoles £ 
, for which there is a dense sampling of the u-v plane. For three angular multipoles 

£ — 1405, 1602 and 1876, this sidelobe suppression along with a low order polynomial 
fitting completely results in residuals of (< 0.02 mK 2 ), consistent with the noise at the 
3cr level. Since the polynomial fitting is done after estimation of the power spectrum it 
can be ensured that the estimation of the HI signal is not biased. The corresponding 
| 99% upper limit on the HI signal is xmb < 2.9, where Shi is the mean neutral fraction 

and b is the bias. 
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1 INTRODUCTION 2.326 GHz (Haslam et al. 1982; Reich 1982; Reich & Reich 

r , , .„ , „„ ,. . , , , 1988; Jonas, Baart, & Nicolson 1998) show that the GDSE 

Observations ol rcdshiitcd 21 cm radiation from neutral hy- , , . , , „ „ ~, . . , 

, TTT . .... . . . . . , has a steep spectral index to be a « —2.8. lhis is also m 

drogen (HI) hold the potential ol tracing the large scale . . " ... .. , , . „, , . 

' , , , broad agreement with the results presented m Platama et 

structure ol the Universe over a large redshift range 20 > , / lnn0 N n . „ , , ,„„„„, , , cnl 

. . , ~ v — al. (1998). Di Matteo et al. (2002) have used the 6C sur- 

2 > 0). This signal (~ mK), is however, buried m the emis- , TT . „ . . . „ ... .„.., , „„„ 

— ' , . , , , , ve y (Hales, Baldwin & Warner 1988), and the 3CR survey 

sion Irom other astrophysical sources which are collectively 

referred to as foregrounds. These foregrounds are dominated 



and the 3 CRR catalogue ( Laing , Riley & Longair 1983) 

to estimate the resolved extragalactic radio sources (point 
by the extragalactic radio sources, the diffuse synchrotron , . .. , 1K „n» TT r, ,. , , 

. . /■/ttvo-o's j n sources) contribution at 150 MHz. Recently Bernardi et al. 



radiation from our own Galaxy (GDSE) and a smaller con 
tribution comes from galactic free-free emission (Shaver et 
al. 1999). The extragalactic point sources have a typical 



(2009) have characterized the power spectrum of the total 

diffuse radiation at 150 MHz at the angular scales of our 

interest. Separating the redshifted HI signal from the fore- 
spectral index ol p ~ —0.8 (Subrahmanyan 2002), with evi- , , . , , , „ , , , 

v , ' grounds, which are several order or magnitude larger (e.g. 

dence of flattening at lower frequencies (Cohen et al 2004) ^ ^ Di Matteo et al. 2002), is currently the 

The analysis ol radio surveys at 408 MHz, 1.42 GHz, and , . , , ,, r m , i r „ i • i i i -r, 

biggest challenge lor 21-cm tomography ol the high redshilt 

Universe. The problem, in principle, can be solved using the 
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t Email:somnath@phy.iitkgp.ernet.in quencies v and v + Av are expected to be uncorrelated 

f Email:saiyad@phys.jdvu.ac.in at separations Av > 0.5 MHz at the angular scales of our 

§ Email:chengalu@ncra.tifr.res.in 



2 A. Ghosh, S. Bharadwaj, S. S. Ali and J. N. Chengalur 



interest whereas the foregrounds, which arise from contin- 
uum sources, are expected to remain correlated over con- 
siderably large frequency separations (Bharadwaj & Sethi 
2001; Di Matteo et al. 2002; Oh and Mack 2003; Di Mat- 
teo et al. 2004; Zaldarriaga, Furlanetto & Hernquist 2004; 
Morales & Hewitt 2004; Bharadwaj and Ali 2005; Santos 
et al. 2005; Wang et al. 2006). We note that the prospects 
of detecting the redshifted 21-cm HI signal are considerably 
better at higher frequencies (e.g., 610 MHz) which probe the 
post-reionization era (z < 6) in comparison to the lower fre- 
quencies (e.g., 150 MHz ) which probe the reionization era. 
Further, the problem of man made radio frequency inter- 
ference (hereafter RFI) is considerably less severe at higher 
frequencies. 

In this paper we report a substantial improvement in 
foreground removal in the context of our earlier work (Ghosh 
et al. 2011, hereafter Paper I) which presents Giant Metre- 
wave Radio Telescope (GMRT 1 Swamp et al. 1991) 610 MHz 
observations towards detecting the post-reionization 21-cm 
signal from z = 1.32. In Paper I we have determined, possi- 
bly for the first time, the statistical properties of the back- 
ground radiation over the angular scales 20 to 10 and 
a frequency band of 7.5 MHz centered at 616.25 MHz. The 
analysis was carried out using the multi-frequency angular 
power spectrum (MAPS) Ct(Av) (Datta, Roy Choudhury 
& Bharadwaj 2007) which jointly characterizes the angular 
I and frequency Av dependence of the fluctuations in the 
background radiation. The measured Ce(Av), which ranges 
from 7 mK 2 to 18 mK 2 , is dominated by foregrounds, the 
expected HI signal being several orders of magnitude smaller 
(Cf l (Av) ~ 10~ 6 - 10~ 7 mK 2 ). The measured signal, for 
a fixed I, is expected to vary smoothly with Av and re- 
main nearly constant over the observational bandwidth. We 
find instead that in addition to a component that exhibits 
a smooth Av dependence, the measured Ct(Av) also has a 
component that oscillates as a function of Av. The ampli- 
tude of the oscillating component is around 1 — 4% of the 
smooth component, and the amplitude and period of oscil- 
lation both decreases with increasing I. We note that similar 
oscillations, with considerably larger amplitudes, have been 
reported in GMRT observations at 153 MHz (Ali, Bharad- 
waj & Chengalur 2008) which is relevant for the signal from 
the reionization era. The origin of the oscillatory signal was 
unclear. 

The oscillatory patterns pose a serious obstacle for fore- 
ground removal (Paper I). It is thus important to identify 
the cause (or causes) and implement techniques to mitigate 
the oscillatory patterns. The fact that the primary beam 
(hereafter PB) pattern changes with frequency across the 
observational bandwidth, not included in our previous anal- 
ysis (Paper I), could be a possible cause. In particular, the 
angular position of the nulls and the side- lobes changes with 
frequency, and a bright continuum source located near the 
null or located in the sidelobes will be seen as oscillations 
along the frequency axis in the measured visibilities. It is 
thus quite plausible that bright sources located near the null 
or the sidelobes of the PB produce the oscillatory pattern in 
the measured Ce(Av) which is estimated from correlations 
amongst the visibilities. One can, in principle, design anten- 
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nas with a frequency-dependent collecting area to produce 
a nearly constant PB pattern. However, to our knowledge, 
none of the upcoming arrays have this feature and we expect 
this issue to be relevant not only for the GMRT but for all 
the arrays planned in the near future. 

The problem can be mitigated by tapering the array's 
sky response with a frequency independent window function 
W(6) that falls off before the first null of the PB pattern and 
thereby suppresses the sidelobe response. It is simplest to 
implement this by multiplying the sky image 1(8) with W(9) 
and using a Fourier transform of this to recalculate the vis- 
ibilities. This, however, will introduce correlations between 
the noise in the recalculated visibilities which is a nuisance 
for estimating Ce(Av). The other option is to work entirely 
with the visibilities in the u — v plane, avoiding the need for 
an image. In this approach the sky response is tapered by 
convolving the visibilities with W(U) - the Fourier transform 
of W{6). As we show later in this paper, it is possible to im- 
plement the convolution without introducing a noise bias in 
the estimated Ce(Av). We note that convolving the visibili- 
ties with a "gridding convolution function" has long been a 
standard practice while making images from interferometric 
data (Section 3., Sramek & Schwab 1989). This convolution 
is done to avoid aliasing which would otherwise occur when 
one uses an FFT to make the image. As discussed in de- 
tail by Sramek & Schwab (1989) the convolution function 
is generally chosen to provide an optimum balance between 
alias rejection and ease of computation. The focus here is 
somewhat different, viz. to use the convolution to strongly 
attenuate the frequency dependent response to the sidelobes 
of the primary antenna pattern. 

RFI sources, which are mostly located on the ground, 
are picked up through the sidelobes. Suppressing the side- 
lobe response is also expected to mitigate the RFI contribu- 
tion. 

Paper I contains a detailed description of the data, 
we mention a few salient features here. The data is taken 
from an archival 30 hours observation centered on 0:2000 = 
12 ,l 36 m 49 s , 5 2 ooo = 62° 17 57" which is situated near Hub- 
ble Deep Field North (HDF-N) . For the present work we 
have analyzed the frequency range 612.5 MHz to 620.0 MHz 
with channels 125 kHz wide. Visibilities were recorded for 
two orthogonal circular polarizations with 16s integration 
time. Calibration was carried out using standard AIPS tasks, 
and the visibilities from the two polarizations were combined 
(V = [Vrr + Vll]/2) for the rest of the analysis. 



2 SIDELOBE SUPPRESSION 

The observed visibilities V = a Cg> I record the Fourier trans- 
form of the sky brightness / convolved with the antenna 
aperture a. The frequency dependence and the sidelobes 
come in through a(U, v) whose Fourier transform gives the 
PB pattern A(6, v). Close to the phase centre, the PB is 
reasonably well modeled by a Gaussian A(9, v) = e~ s 
where the parameter 9o is related to the FWHM of the PB 
as 80 « 0.6 x 0fwhm, and do — 25 .8 (#fwhm = 43 ) at 610 
MHz for the GMRT with 6 oc v' 1 . Note that the Gaussian 
model for A(8, v) breaks down away from the phase center 
where we have the sidelobes and nulls. 



GMRT observations : Improved 21-cm foreground removal 3 



We mitigate the effect of the frequency dependent side- 
lobe pattern by convolving the observed visibilities with a 
suitably chosen function W(U). The sky response of the con- 
volved visibilities V c = W ® V is modulated by the window 
function W(6) which is the Fourier transform of W(U). We 
have used a window function W(8) = e~ e / e ™ with 8 W < 80 
to taper the sky response so that it falls off well before the 
sidelobes. Note that 8 W and W(8) are both frequency inde- 
pendent. We parametrize 9 W as 9 W — f#o with f < 1 where 
80 here refers to the value at the fixed frequency 610 MHz. 

We can evaluate the convolved visibilities V c on a grid 
in u - v space (V c (Ui,v) = W(U l - U a )V(U„,i/)) 
and use these to determine the two-visibility correlation 
defined as V 2 (Ui,Av) = V c (Ui, i>)V*(Ui, v + Av). Here 
U = (it, v) refers to a two-dimensional baseline, Uj refers 
to points on the grid in u — v space and U a refers to 
the different baselines in the observational data. This way 
of estimating 14(Uj, Av), however, introduces a positive 
noise bias (e.g. Begum et al. (2006)) which is not desir- 
able. We use, instead, the estimator V2(Ui, Av) = K^ 1 x 
[W(Ui - U„) W*(U t - U 6 )V(U a , v)V*(U b , v + Ai/)] , 

where K = £ a _^ W(Ui - U a )W'*(U l - U b ) is a normal- 
ization constant. The noise bias is avoided by dropping 
the self-correlations (i.e.. the terms with a = b). We have 
used a grid of spacing AU g = Vla2/(Tv8 w ) = 0.265#^ 1 
which corresponds to half of the FWHM of W(U), and 
we have estimated V^(Ui,Ai/) at every grid point using 
all the baselines within a disk of radius 2 AU g centered on 
that grid point. We finally determine Ct(Av) using (Ali, 
Bharadwaj & Chengalur 2008) 



V 2 (U,Av) = ^ 'c e (Av)Q(Av). 



(1) 



where £ = 2irU, 8' 2 = 8 2 + 6 W 2 and Q{Av) is a slowly 
varying function of Av which accounts for the fact that we 
have treated 9' 2 and (8I v /dT) as constants in our analy- 
sis. We have used Q(Au) = 1 which introduces an extra 
Av dependence in the estimated Ce(Av) . This, we assume, 
will be a small effect and can be accounted for during fore- 
ground removal. The data has been binned assuming that 
the statistical properties of the signal are isotropic in U. 



2.1 Simulation 

In order to demonstrate the efficacy of our technique of side- 
lobe suppression using a simulated data set, the GMRT an- 
tenna was modeled as a circular aperture of D = 45 m diam- 
eter with a circular disk of diameter Di = lmat the center 
blocked due to the feed. This gives a normalized antenna 
beam pattern 



A v {8) = 



(D 2 -D2)2 



D 



■ Ji (tt8D/\) 
(tt8D/X) 



3 Ji(7rflPi/A) 
1 (tv8D 1 /X) 



(2) 



where Ji is the Bessel function of the first kind of order one. 
Several point sources were randomly placed in the region 
near fifth null of this beam pattern. The expected number 
of point sources and their flux distribution can be predicted 
from the differential source counts (Garn, Green, Riley et al. 
2008) . Given that the present simulation has the limited aim 
of demonstrating the efficiency of our sidelobe suppression 
technique, we have instead used a cartoon model with five 



hundred point sources each having the same flux density 
150 mjy. The value of the flux density was chosen so that the 
simulated Ct(Av) has a value around ~ 10 mK 2 , comparable 
to our measured Ce(Av) (Figure 7, Paper I). 

Next, we have randomly generated baselines, and cal- 
culated the corresponding visibilities for the five-hundred 
point sources. It may be noted that the simulated baselines 
have a density of 0.3/m 2 which is roughly consistent with 
the GMRT central square. The antenna beam pattern (eq. 
2) was used in calculating the simulated visibilities. We have 
used the simulated visibilities to estimate Ce(Av) . Then, we 
have fitted a third order polynomial in Av to subtract out 
the component of the sky signal that varies slowly with fre- 
quency. The residual Ct(Av) has an oscillatory pattern with 
amplitude (~ 10~ 2 mK 2 ). The procedure was repeated after 
introducing a convolution with / = 0.8. This tapers the sky 
response to ~ 0.8% of its peak value at the first null of the 
beam and considerably reduces the response to the sources 
beyond the first null. We find that the residuals, after poly- 
nomial subtraction, are in the range 10 -5 mK 2 to 10~ 4 mK 2 , 
i.e. the residuals are suppressed by a factor of around 10 2 
to 10 3 . Further, after convolution the residuals appear to be 
noise like and do not show any noticeable oscillatory feature. 
Note that although we have shown results of our simulation 
for only a single I value, the results are very similar at other 
i values. 



3 RESULTS AND CONCLUSIONS 

We first investigate whether sidelobe suppression at all helps 
to reduce the oscillations that were reported in Paper I. This 
is quantified using the dimensionless decorrelation function 
Kt(Av) = Ce(Av)/Ce(0). We have considered f = 0.4,0.65 
and 0.8 which respectively correspond to a tapered sky re- 
sponse with FWHM 17'. 2, 28' .0 and 34'. 4 as compared to 
the GMRT PB which has a FWHM of 43' at 610 MHz. The 
results are shown in Figure 2 for the four smallest I values 
for which we have measured the binned Ce(Av) . The oscil- 
latory patterns are distinctly visible in the cases where the 
tapering has not been applied. We see that in most cases the 
oscillations are considerably reduced and are nearly absent 
when tapering is applied. We do not, however, notice any 
particular qualitative trend with varying f . The oscillatory 
pattern with a large period in Av seen at I = 1405, the 
smallest £ value, persists even after the tapering is applied. 
This, however, does not pose a problem for foreground re- 
moval as the Ce(Av), after tapering, is well fitted by a low 
order polynomial and can be successfully removed. 

The tapering, which has been implemented through a 
convolution, is expected to be most effective in a situation 
where the u — v space is densely sampled by the baseline 
distribution. Our results arc limited by the patchy u — v 
coverage of the present observational data. This also pos- 
sibly explains why the oscillations persist after tapering at 
the smallest I value. Tapering the the field of view has a 
drawback in that this increases the cosmic variance which 
scales as f _1 . It is possible to compensate for this by repeat- 
ing the entire analysis after adding phases to the visibilities 
so as to shift the center of the field of view. We have not 
attempted this here, and our entire analysis is restricted to 
a single field of view. The increase in the cosmic variance 
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Ai/ MHz Ai/ MHz 

Figure 1. The data points in the top panels show C((Au) estimated without (left) and with (right) sidelobe suppression. The solid lines 
show the respective best fit polynomials used for foreground subtraction. The residuals, after foreground subtraction, are shown in the 
corresponding bottom panels. 



can also be used as a guiding principle in choosing the value 
of f. It is most advantageous to use the smallest value of 
f where the oscillations are adequately suppressed. Of the 
three values of f that we have considered here, we find that 
foreground subtraction is most effective for f = 0.65, and we 
use this for the entire subsequent analysis. 

The measured Ci(Av) shown in Figure 3 is foreground 
dominated and it is seen to vary smoothly with increas- 
ing Av. On the contrary, the predicted contribution from 
the HI signal C™ (Av) decreases very rapidly with increas- 
ing Av (Bharadwaj & Sethi 2001). We have used a polyno- 
mial fitting technique (Paper I) to identify and subtract out 
any smoothly varying component from the measured Ci(Av) 
and thereby remove the foreground contribution. The possi- 
bility that along with the foregrounds the fitting procedure 
may also remove a part of the signal is a major concern. For 
each £ value we have estimated Avo.i which corresponds to 
the frequency separation where C™(Av) first falls to less 
than 10 % of the peak value C™(0). The bulk of the HI sig- 
nal is localized within Av < Avo.i which is excluded when 
estimating the slowly varying foreground contribution. We 
have used the range Avo.i < Av < 6 x Avo.i to estimate 
the coefficients of the best fit 4 th order polynomial and we 
use this to subtract out the foreground contribution from the 
entire range Av < 6 x Avo.i- The value of A^o.i and the 
polynomial fit are both shown in Figure 3. Note that Afo.i 
decreases with increasing I, and it is less than 0.5 MHz for 
all the I values that we have considered. Tests with simula- 
tions (Paper I) show that very little of the HI signal is lost 
in this subtraction procedure. 

It is noteworthy that a variety of foreground removal 
techniques (e.g. Mcquinn et al. 2006; Jelic et al. 2008; Gleser 
et al. 2008; Harker et al. 2009; Liu, Tegmark & Zaldarriaga 
2009; Liu et al. 2009; Petrovic & Oh 2011) all attempt to 
remove the foregrounds from images or visibilities before de- 
termining the power spectrum. The HI signal is spread out 
along the entire frequency axis and all these techniques run 
the risk of removing a part of the HI signal along with the 
foregrounds. On the contrary, we have attempted foreground 



removal after determining the angular power spectrum. The 
signal, here, is localized in Av , and it is possible to re- 
duce the risk of removing a part of the HI signal by suitably 
tuning the fitting procedure. 

We expect the residual Cf' ES (Av) that remains af- 
ter polynomial subtraction (Figure 4) to contain only the 
HI signal and noise provided the foregrounds have been 
successfully removed. In our observations the HI signal is 
much smaller than the noise a ~ 0.01 mK 2 . We thus ex- 
pect the residuals to be consistent with noise provided the 
foregrounds have been completely removed. We find that 
the residuals are consistent with noise at the 3 a level (ie. 
Cf ES (Av) < ± 3 it) at the three smallest values of £ 
(1405, 1602, 1876). This establishes that our foreground re- 
moval technique works. This technique, however, is not as 
successful at £ = 2214 where a single point of Av — has a 
value that is somewhat larger than 0±3 a. We have also car- 
ried out the entire analysis for several other, larger, values 
of £ for which the results have not been shown here. While 
sidelobe suppression works quite well in removing the oscil- 
lations, we are unable to completely remove the foregrounds 
for the other I values. 

The fact that we have three £ values where the fore- 
grounds have been completely removed and Cf ES (Av) is 
consistent with noise allows us to place an upper limit on 
the HI signal. We assume that the HI traces the dark matter 
with a possible linear bias b whereby C™(Av) can be calcu- 
lated (Paper I) in terms of the dark matter power spectrum, 
and the cosmological parameters for which we have used the 
values (fi m0 , Qao, h,a 8 ,n s ) = (0.3, 0.7, 0.7, 1.0, 1.0). The HI 
signal is now completely determined upto a proportionality 
factor C™(Av) oc [xhi&] 2 where xm is the mean neutral 
fraction of hydrogen gas. We have considered xmb as a free 
parameter, and performed a likelihood analysis using our 
observational data to place an upper limit on xmb. In our 
analysis we have assumed that the Ce(Av) at the differ- 
ent I values are independent. For each £, we have only used 
Ce{Av) in the range Av < Avo.i; the Av values that were 
used to estimate the polynomial for foreground removal were 
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Figure 2. The measured ki(Au) as a function of Au for the £ value as shown in each panel. The results from our earlier analysis (Paper 
I) and for a tapering f = 0.4, 0.65, 0.8 arc shown by the thick solid, dashed, thin solid and dotted curves respectively. 
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Figure 3. The data points show the measured Ci(Au) as a function of Au for the I value shown in each panel. The error-bars indicate 
the 3<r system noise. The solid curve shows the 4 th order polynomial fits and the dotted vertical line shows Auq.i- 



excluded in the likelihood analysis. The covariance between 
Ci(Au) at different Au values was taken into account in es- 
timating the likelihood. We place an upper limit xmb < 2.9 
with 99% confidence using our observational data. In Paper I 
a high-pass filter had been applied to remove the oscillatory 
pattern at a single I value (the smallest) to give an upper 
limit of 7.95 for xmb. Table 1 summarizes the results of the 
current work, and presents a comparison with the results of 
Paper I. It is clear that the present analysis is a consider- 
able improvement over the previous result. In conclusion we 
show that sidelobe suppression can bring about a consider- 
able improvement in foreground removal. 
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Table 1. We show ACt {Au = 0), the 1 — c error, for the I values 
where it was possible to successfully remove the foregrounds. The 
upper limits on xmb from the single I value in Paper I and from 
combining the three different t values in the current work are also 
shown. 
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Figure 4. The residuals are shown with 3<r error bars (system noise only). 
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